 function band = nwlag_SPJ(y , x , alpha , omega, kernel)
 %This code chooses the NW lag following Sun et al. (2008 Econometrica)
 % y -- dependent variable
 % x -- independent variable
 % alpha -- two-sided critical value
 % omega -- relative weight on type I errors
 % kernel -- 1: Parzen kernel; 2: QS kernel
 %% Check dimensions
 [T , ~] = size(x);
 
 if size(y , 1) ~= T
     error('Wrong dimension')
 end
 
 %% Values
 if kernel == 1 %Parzen
   g = 6;
   c = 0.539;
 elseif kernel == 2 %QS
   g = 1.421;
   c = 1;
 else 
     error('Kernal type wrong')
 end

 z_alpha_s = (norminv(1 - alpha/2))^2;
 
 %% First estimate the data with OLS
 b_x = olsgmm(y , [ones(T , 1) , x] , 0 , -1);
 u = y - [ones(T , 1) , x] * b_x;
 
 %% Fit AR(1) to residuals
 b_rho = olsgmm(u(2 : end) , [ones(T - 1 , 1) , u(1 : end - 1)] , 0 , -1);
 
 %% Paramters
 d_hat = 2 * b_rho(2)/(1 - b_rho(2))^2;
 delta = 2;
 
 %% Calculate quantities
 Q = d_hat * (omega * ncx2pdf(z_alpha_s , 1 , 0) - ncx2pdf(z_alpha_s , 1 , delta^2));

 J = 1000;%Big enough
 K_delta_j = nan(J , 1);
 z_K = z_alpha_s;
 for j = 1 : J
     K_delta_j(j) = (delta^2/2)^j * exp(-delta^2/2) * (1/cumprod(j)) * j/(gamma(j + 0.5) * 2^(j + 0.5)) ...
                  * z_K^(j - 0.5) * exp(-z_K/2)/z_K;
 end
 
 K_delta = nansum(K_delta_j);
 %% Determine bandwidth
 if Q > 0
     band = (2 * g * Q/(c * z_alpha_s * K_delta))^(1/3) * T^(-2/3) * T;
 else
     band = log(T);
 end
 end